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ABSTRACT 

This  thesis  addresses  the  problem  of  identifying  the 
dynamics  of  a  linear  system  in  the  frequency  domain.  An 
alogrithm  operating  on  the  Fast  Fourier  Transform  (FFT)  of 
blocks  of  signals  is  developed  and  its  performance  evaluated 
through  computer  simulations.  Several  properties  are  tested, 
in  particular,  its  convergence  and  its  capabilities  of 
identifying  the  frequency  response  of  the  unknown  system. 
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I .  INTRODUCTION 

The  problem  of  identifying  the  dynamics  of  a  physical 
system  from  its  input-output  behavior  is  a  fundamental  issue 
in  modeling  and  control.  The  most  general  approach  is  based 
on  the  assumption  of  a  mathematical  model,  often  linear  and 
time  invariant,  and  its  parameters  are  determined  by 
optimization  of  an  error  criterion. 

A  wide  class  of  identification  algorithms  exists  in  the 
literature.  Several  books  have  been  written  on  this  matter, 
see  for  example  [Ref.  1,  2].  Of  particular  interest  are  the 
recursive  algorithms,  where  the  estimated  parameters  are 
updated  on-line  based  on  new  input-output  data  collected  at 
each  sampling  instant.  In  most  of  the  approaches  on  recursive 
identification  the  parameter  estimates  are  updated  for  every 
new  set  of  data  points;  as  a  consequence,  the  estimate  at  any 
time  is  updated  up  to  the  current  data. 

A  different  approach  is  based  on  Block  Processing  (BP)  , 
where  the  parameters  are  recursively  updated  on  the  basis  of 
the  data  within  nonoverlapping  intervals  of  time.  Although 
this  introduces  a  delay  between  the  time  we  measure  the  input- 
output  data  and  the  time  we  use  it  in  the  estimation,  it  has 
the  advantage  that  time  averaging  within  the  time  interval 
reduces  the  effect  of  disturbances. 


Block  Processing  techniques  have  been  investigated  by 
several  authors.  Of  particular  interest  for  this  thesis  is 
the  approach  developed  by  Mansour  [Ref.  3]  in  the  adaptive 
filtering  contest,  where  the  parameters  are  updated  on  the 
basis  of  the  DFT  (Discrete  Fourier  Transform)  of  the  data 
within  each  time  interval. 

In  this  thesis  an  adaptive  identification  algorithm 
operating  on  the  DFT  (Discrete  Fourier  Transform)  of  blocks  of 
input  and  output  data  is  introduced.  In  particular,  the 
algorithm  identifies  the  frequency  response  of  the  system's 
linear  model  by  applying  several  recursive  estimation 
techniques  and  adapting  them  to  the  frequency  domain  approach. 
We  will  investigate  estimation  methods  based  on  the  Projection 
Algorithm  and  Recursive  Least  Squares  [Ref.  4]  and  the  results 
will  be  compared.  Even  though  Recursive  Least  Squares  is  a 
complex  and  time  consuming  method,  it  is  preferable  when 
accuracy  and  speed  of  convergence  are  needed. 

This  thesis  is  organized  as  follows.  Chapter  II 
introduces  the  concept  of  Block  Processing  and  develops  two 
algorithms  for  recursive  identification  in  the  frequency 
domain.  Proof  of  convergence  is  also  part  of  this  chapter. 
Simulation  studies  are  introduced  in  Chapter  III  showing  the 
effectiveness  of  the  algorithms  and  the  conclusion  is  in 
Chapter  IV. 


II.   IDENTIFICATION  OF  LINEAR  SYSTEMS 


A.   STATEMENT  OF  THE  PROBLEM 

Consider  an  LTI  (Linear  Time  Invariant)  system  with  input, 
output  signals  given  by  u(t)  and  y(t)  respectively.  The 
problem  we  address  is  the  estimation  of  the  parameters  of  a 
linear  model  in  order  to  predict  future  values  of  the  output 
y  given  past  measurements  of  the  input  and  output  signals. 
The  general  block  diagram  is  given  in  figure  2.1.  The  block 
z"1  denotes  time  delay  by  one  clock  pulse. 
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Figure  2.1  Block  diagram  of  the  system 


This  figure  shows  the  general  scheme  where  y(t)  is  the 
prediction  of  y(t)  based  on  the  past  input,  output,  and  the 
estimated  model.  The  purpose  is  to  determine  a  model  for  the 
plant  so  that  the  error  e(t)  is  miniminized.   In  particular, 


the  model  is  restricted  to  be  an  LTI  system  in  either  one  of 

the  two  forms: 

ARMA  (Auto  Regressive  Moving  Average) 

y(t)  =a1y(fc-l)  +.  .  .+any(t-n)  +jb1u(fc-l)  +  .  .  .bnu(t-n)     (2.1) 

where  ai ,  b.  are  constants;  or 

CONVOLUTION  form,  where  the  output  input  signals  are  related 

by  the  convolution  sum. 

y(t)=h(t)  *u(t)  (2.2) 


y{t)=Y,  h(x)uit-x)  (2.3) 


The  problem  of  adaptive  freguency  response  identification 
mainly  has  been  concentrated  on  the  approximation  of  infinite 
impulse  response  systems  by  finite  impulse  response  models. 
In  this  thesis,  we  will  consider  the  identification  of  stable 
plants  only,  for  which  the  impulse  response  h(t)  decays  to 
zero  as  time  increases.  For  this  reason,  the  infinite  impulse 
response  h(t)  can  be  approximated  by  a  finite  seguence  with 
length  N.  In  this  way,  the  moving  average  (MA)  model  of  the 
following  form  is  obtained. 


y(t)  =h(0)  u(t)  +h(l)  u(t-l)  +-+h(N-l)  u(t-(N-l)  )    (2.4) 

where  h(t),  for  0  <  t  <  N-l,  is  the  truncated  impulse 
response,  We  choose  the  block  length  N  to  be  a  power  of  2  in 
order  to  apply  FFT  (Fast  Fourier  Transform)  techniques. 

The  approach  we  consider  is  based  on  a  frequency  response 
framework.  In  particular,  the  estimated  model  is  computed 
from  the  Fourier  Transform  of  blocks  of  data  collected  from 
measurements  of  the  input  and  output. 

There  are  several  reasons  which  motivate  this  approach. 
In  many  cases  of  interest,  a  model  which  is  accurate  on  a 
desired  frequency  range  only  is  identified  for  the  plant.  For 
example,  if  we  want  to  model  the  low  frequency  spectrum  of  the 
system  with  the  frequency  response  approach,  we  select  weights 
in  order  to  select  frequencies  of  interest. 

In  the  next  section,  a  recursive  algorithm  for  the 
identif icaton  of  the  frequency  response  is  introduced. 

B.   RECURSIVE  FREQUENCY  DOMAIN  APPROACH 

Block  Processing  (BP)  techniques  for  on  line 
identification  of  linear  models  have  been  investigated  by 
several  authors  [Ref.  3].  With  this  approach  the  estimated 
model  parameters  are  sequentially  updated  on  the  basis  of 
blocks  of  data,  rather  than  at  each  data  point.  In  this 
section,  we  investigate  a  BP  adaptive  algorithm  which  operates 


on  the  Fast  Fourier  Transform  (FFT)  of  blocks  of  input  and 
output  data. 
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Figure  2.2  Blocks  of  Input  and  Output  Data 

There  are  some  advantages  in  dealing  with  data  by  blocks 
rather  than  individually.  Block  processing  averages  the  data 
within  the  block,  resulting  in  better  noise  rejecton.  In 
order  to  illustrate  the  BP  approach  in  the  ARMA  model  case, 
consider  the  model  in  (2.1)  and  the  data  yk(t),  uk(t)  within 
the  k-th  time  block,  i.e. 


yk(t)  =y(kN+t) 


uk(t)  =u(kN+t)  (2.5) 

for  t  =  0,  1,  ...  ,  N-l.   If  the  block  size  N  is  larger  than 
the  system  order  n,  we  can  write  (2.1)  as 

yk(t)  =a1yk(t-l)  +a2yk(t-2)  +b1uk(t-l)  +b2uk(t-2)         (2.6) 


for  a  second  order  case  with  2  <  t  <  N-l.  It  is  easy  to  see 
that  (2.6)  can  be  written  in  convolution  form  as 

yk(t)  =a(t)®yk{t)  +b(t)®uk(  t)         nztzN-1  (2.7) 

with  a  and  b  sequences  in  Rn 

a=[0   ax  a2  -  an  0   •••  0] 

b=[0  bx  b2  •••  bn  0   -   0]  (2.8) 

and  ®  denoting  circular  convolution. 

In  the  approach,  the  well  known  relationship  between  the 
DFT  and  the  circular  convolution  in  the  frequency  domain  is 
exploited.  By  defining  {  Yk(l)  ,  1=0,1,  ...,N-1  }  and  {  Uk(l) 
,  1=0,1,  ...,N-1  }  as  the  DFT  of  the  output  and  input  data 
yk(t)  ,uk(t)  in  the  k-th  block,  we  can  define   Y  k(l)  as 

Yk{l)  =A{1)  YA1)  +B{1)  UA1)  2=0,1,-,  N-l.  (2.9) 


Clearly  from  (2.7)  and  (2.9),  we  can  write  the  important 
relationship 

yk(t)=yk(t)  nztiN-1  (2.10) 

with   y  k(t)  =  IDFT[  Y  iJ1)]'  and  IDFT  denoting  inverse  DFT. 
A  linear  estimation  algorithm  for  A(l)  and  B(l)  in  (2.9) 
can  be  determined  by  defining  an  error  signal 


ek(t)  = 


yk{t)-yk(t)  nztzN-1  (2    11} 

0       0^t^72-l 


where  y  k(t)  =  IDFT[  A  k(l)Yk(l)  +  Bk(l)Uk(l)  ,  1=0,1, ...  ,N-1  ]  and 
Ak,  B  k  are  estimates  available  at  the  k-th  block.  From  the 
above  definitions,  we  can  write  the  recursive  estimates  for  A, 
B  by  the  algorithm.   Let 

$k(l)=[Ak(l),BkU)]T 

Xk(l)=[Yk(l)  ,Uk(l)]T  (2.12) 

for  1=0,  1,  ...,  N-l,  k  =  0,  1,  2,  ...;  also  let 

Ek(l)  =DFJ[ek(t)]  (2.13) 

with  e.  as  in  (2.11).   Then  the  recursion  is  such  that 


v  ,is    A  ,11       X'k(l)Ek(l) 
^1(J)=©jc(J)+__ £ _       (2.14) 

i+m5x^(i)i2+i^d)i2) 


For  this  algorithm  we  can  show  the  following: 

(a)  ||  6k+1(l)  -  6(1)  ||2  <  ||  0k(l)  -  0(1)  ||2  (2.15) 

for  all  1=0,  1,  2,  .  ..,  N-l,  and  for  all  k; 

(b)  If  the  input  and  output  data  uk(t)  ,  yk(t)  are  bounded  for 
all  k  and  t,  then 

£i~  ek(t)=0  for  all  OztzN-1  (2.16) 


In  this  result  we  see  that  the  parameter  error  decreases 
with  time  (eqn.  (2.14))  and  the  error  between  the  model  and 
the  plant  ek(t)  tends  to  zero,  i.e.  the  model  output  tends  to 
follow  the  plant  output. 

Proof:  Let  us  define  the  parameter  error 

8k(l)=Qk(l)-Qk(l)  (2.17) 

Then  from  (2.14)  and  (2.17)  we  can  write  the  recursion 

&k^(l)  =&k(l)  -\ik(Xk(l)  Ek(l)  )  (2.18) 


where 


^  l+m**(pk(l)f+\YkU)\2)  (2*19) 


By  taking  the  magnitude  sguare  of  both  sides  of  (2.18),  we 
obtain 

=  [©^(2)  -^^(DiTjd)]  [©^(J) -11^(2)^(2)] 

=||©/t(2)||2-nJ^;r(2)^r(2)©^(2)-jiJce;T(2)xJ;(2)4(2) 

+  \i2kE*kTx£(l)X*k(l)Ek(l)  (2.20) 


By  taking  the  sum  of  all  frequency  components  in  both  sides  of 
(2.20),  we  obtain  the  expression 


EP*»u>l,-Efc«>r-i»iE 

i=o  i=0  J«0 


El^i(^l2=EI^(J)|2-^E  B?ii)xlu)**u) 


1=0 


N-l 
2=0 


+  [!*£  »&*{!)  XZU)  Xk(l)  4(2)  (2.21) 


Applying  Parseval ' s  theorem,   we  can  relate  data  in  the 
frequency  domain  to  the  corresponding  time  domain  as 
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£  [E*/{1)]   [xZ(l)&k(l)]=y£§k(t)tk(t)  (2.22) 

j=o  t-o 


N-l  N-l 

^ek(t)tk(t)=Y: 

C=0  C=2 


E**(t)eJt(t)-EpJt(C)|a^O  (2.23) 


where  we  define 


eIr(t)=J£>FTUJ(J)©t(I)]  (2.24) 


.^  v  ui  -  j.ur  i  is\k   v  j.  )  ^>k 

Equation  (2.23)  comes  from  the  fact  that 

§k(t)=0  Oztzn-1  (2.25) 

ek(t)=yk(t)  -yk(t)  =ek(t)  nztzN-1  (2.26) 

and 

ek(t)=ek(t)  n<t<N-l  (2.27) 

Still  using  Parseval's  theorem,  we  can  write 

ElP*u>IW  (2-28) 


ii 


After  some  manipulations  and  the  fact  that  /uk|Xk(l)|2  <  1  for 
all  k  and  1  =  0,  1,  ...,  N-l,  we  can  bound  the  parameter  error 
at  the  end  of  the  (k+l)th  block  as 


'*♦! 


r  *  w2-2m^mw  (2-29> 


and  therefore 

P*.J2  *  Pf-iHfit  <2-30> 

Finally,  since  ||  ©k  II  ^  0  for  all  k  and  is  a  nonincreasing 
sequence  as  in  (2.30)  ,  the  increment  /liJ|  e  k||2  must  tend  to  zero 
as  k  tends  to  infinity.   Therefore 

lim llg*f =0  (2    31) 

*"°°  l+MAJ§Yk(l)\2+\Uk(l)\2} 

which  proves  the  result. 

The  significance  of  this  result  is  that  we  can  estimate 
the  parameters  of  the  given  system  (in  terms  of  Ak(l),  Bk(l)) 
based  on  the  frequency  content  of  each  block  of  data. 
Furthermore,  the  estimation  is  recursive,  and  the  convergence 
of  the  prediction  error  e(t)  to  zero  is  guaranteed  at  least  in 
the  ideal  case.  It  will  be  shown  in  the  simulations  that  even 
in  the  presence  of  measurement  noise,  the  convergence  of  the 
error  to  small  values  is  still  satisfactory. 
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In  the  implementation  of  the  estimation  algorithm 
introduced  above,  we  need  to  compute  the  DFT  of  the  "windowed" 
error  term  ek(t)  in  (2.25)  and  (2.26).  We  can  see  that  in  the 
implementation  two  operations  of  IDFT  are  required  to  compute 
yk(t)  .  Using  DFT  again,  Ek(l)  is  obtained,  which  is  used  in 
the  recursion  (2.14). 

An  algorithm,  which  does  not  require  this  sequence  of  IDFT 
and  DFT  in  its  implemention,  is  introduced  next.  The  new 
algorithm  also  offers  the  advantage  of  being  able  to  use 
different  recursive  estimation  techniques  (such  as  the 
Recursive  Least  Squares)  which  exhibits  faster  convergence. 
However,  this  is  obtained  at  the  expense  of  added  complexity. 

In  order  to  introduce  the  technique,  let  us  recall  the 
signal 

yk  (t)=a(  t)  ®yk  ( t)  +Jb  ( t)  ®uk  ( t)  (2.32) 

with  a,  b  being  the  plant  parameters,  and  the  equality 

yk(t)=yk(t)  nztzN-1  (2.33) 

Going  to  vector  notation,  let  us  write  (2.33)  in  the  following 
form 
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yk- 


yk(n) 

yk(N-2) 
yk(N-i) 


[o]  [0] 
[0]  [J] 


yk(o) 

ykin) 

yk(N-2) 
yk(N-i) 


=hyk 


(2.34) 


where  [0]  and  [I]  denote  blocks  of  zeros  and  the  identity 
matrix  respectively.   Definition  of  the  Fourier  matrix  F  as 


Fm,n=e 


• ,   2nmn > 


(2.35) 


allows  the  computation  of  the  DFT  of  a  sequence  as  a  matrix 
operation.   In  this  way,  we  obtain  from  (2.34) 


Fyk=Fhyk 


(2.36) 


and  therefore,   after  simple  manipulations,   the  following 
equation, 


Fy^iFhF'1)  (Fyk) 


(2.37) 


It  is  easy  to  see  that  the  vectors  Fyk  and  F  y  k  are  the  arrays 
of  DFT  coefficients  of  the  respective  sequences  yk  and  yk, 
which  yields 


Yk=HYk 


(2.38) 


Recall  from  (2.32)  that 
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YA1)=AA1)  YA1)  +BADUA1) 


k^J-'  "k 


(2.39) 


for    1=0,     1,     2,     .  ..,    N-l.     Let    us    define    YU    as    the    following 
expression 


YU= 


Yk(0)         0 

0         Yu(l)    ■■■ 


0 


0 


i     Uk(0)  0 

i  0  Uk(l) 

■         :  0 


0    Y.(N-l) 


0 
0 

0    U.iN-1) 


(2.40) 


and  combine  (2.38)  and  (2.39)  to  yield 


AAO) 


Yk=  (HYU) 


Ak(N-l) 

Bk(0) 
BAN-1) 


(2.41) 


The  recursive  estimation  algorithm  is  derived  from  (2.41) .  In 
particular,  we  know  that  (2.41)  can  be  broken  into  a  sequence 
of  linear  equations  of  the  form 


Yk(l)=$l(l)@ 


1=0,1,2, . . . ,N-1 


(2.42) 


where 
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0=[^(O)  ...  Ak(N-l)    Bk(0)     -  Bk(N-l)]T  (2.43) 

is  the  vector  of  unknown  parameters,  and 

*J(2)=[iJ(2,0)  -  H(1,N-1)]  (YU) 

=  [H(l,0)Yk(0)    -  H(1,N-1)  YAN-1) 


,#(2,0)  £7.(0)  -  H{1,N-1)UAN-1)] 


(2.44) 


In  (2.44)  the  coefficient  H(i,j)  are  constant  (in  the  sense 
that  they  do  not  depend  on  the  block  k)  and  are  determined  by 
the  window  matrix  h  only. 

From  (2.42)  we  can  use  any  recursive  algorithm  to  estimate 
0.  In  particular,  a  Recursive  Least  Squares  [Ref.  2] 
algorithm  applied  to  complex  data  is  going  to  yield  the 
following  recursion 


(2.45) 
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1*1    i  p£®k(D®k(i)Pk 
H  -pi-      T       i  (2-46) 


for  1  =  0,  1,  2,  ...,  N-l,  with  initial  conditions  at  the 
beginning  of  each  time  block  given  by 

K=K-i      S      P°k=°2I  (2.47) 

with  I  the  identity  matrix  and  a2  an  arbitrary  constant 
parameter.  In  general,  the  constant  a  is  chosen  to  be  a  large 
value,  in  comparison  with  the  order  of  magnitude  of  the 
parameters  6 . 

Although  far  more  complicated  to  implement,  the  algorithm 
(2.45),  (2.46),  and  (2.47)  has  better  convergence  properties 
than  the  one  shown  previously. 

In  the  next  section  the  two  algorithms  are  applied  to  the 
identification  of  a  linear  time  invariant  system. 
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III.   SIMULATION  STUDIES 

A.  INTRODUCTION 

The  estimation  techniques  discussed  in  the  previous 
chapter  have  been  simulated,  and  used  to  estimate  the 
frequency  response  of  a  system.  In  addition,  two  computer 
programs  written  in  MATLAB  have  been  developed.  The  first  one 
is  the  main  program  given  in  Appendix  A  which  simulates  the 
entire  system  in  an  iterative  fashion.  The  main  program 
estimates  the  value  of  the  frequency  response  of  a  given 
system  at  a  particular  frequency  and  compares  it  with  the 
corresponding  value  of  the  original  system. 

The  second  program  is  a  function  type  of  subroutine  which 
implements  the  Recursive  Least  Squares  algorithm  [Ref.  2]. 

The  proposed  identification  method  has  been  tested  using 
several  different  inputs  in  the  program  in  order  to  simulate 
and  analyze  the  effect  of  various  excitaton  conditions.  In 
the  next  section,  a  first  order  recursive  difference  equation 
is  used  as  an  example. 

B.  SIMPLE  DERIVATION 

The  system  being  used  as  an  example  is  a  first  order 
system  described  by  the  difference  equation. 
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y(fc)  =0.3y(fc-l)  +5.  Ou(fc-l)  (3.1) 

In  like  manner  (3.1)  can  also  be  written  in  terms  of  the 
impulse  response, 

oo 

y(t)  =-h(0)u(0)+h(l)u(t-l)+.  .  .  =  £  h(m)u(t-m)         (3.2) 


Since  the  system  is  stable,  its  impulse  response  decays  to 
zero,  and  (3.2)  can  be  approximated  by  a  finite  sum 

y(fc)=h(0)u(t)+h(l)u(t-l)+. . .+h(N-l)  u(t-N+l)        (3.3) 

In  the  followig  simulations  N  is  chosen  as  N  =  32  data  points. 
The  transfer  function  H(z)  of  the  original  system 

tf(z)=XLEL  = § (3.4) 

U(z)      z-0.3 


yields  its  frequency  response 


H(eJQ)=-— -5 (3.5) 

ej6-0.3 


Due  to  the  use  of  FFT  methods,  we  estimate  the  frequency 
response  at  discrete  frequency  points  6  =  27rl/N,  1  =  0,  1,  2, 
...,  N-l.   For  this  purpose,  let  us  define 
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-j(-S^S) 


Hd)=j2h^n)e      K  (3*6) 

and   its   finite   approximation 

N-l  _.,  2Kln) 

H{l)=Y,h(n)e         N  (3.7) 


n=0 


The  next  step  is  to  calculate  the  frequeccy  response  of  the 
original  system. 

C.   STRUCTURE  OF  THE  PROGRAM 

The  algorithms  introduced  in  the  last  chapter  have  been 
tested  in  the  identification  of  a  discrete  time  system.  In 
particular  we  look  at  the  problem  of  identifying  the  frequency 
response  of  the  system  by  fitting  a  finite  impulse  response 
system  to  the  input-output  data. 

The  software  developed  consists  of  two  programs:  a  main 
program  which  produces  the  input  output  data  used  in 
identification,  and  a  subroutine  which  implements  the 
Recursive  Least  Squares  identification  shown  in  Chapter  II. 

The  effectiveness  of  the  estimation  is  assessed  on  the 
basis  of  two  criteria:  the  error  between  actual  output  of  the 
system  and  predicted  output,  and  the  error  between  the 
frequency  response  of  the  system  and  the  frequency  response  of 
its  estimate. 
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Since  we  try  to  fit  a  finite  impulse  response  model  to  a 
recursive  system  (with  impulse  response  of  infinite  length) , 
we  cannot  pretend  to  estimate  the  frequency  response  at  all 
frequencies  in  an  effective  way.  Therefore,  we  input  test 
signals  at  finite  number  of  frequencies  and  we  look  at  the 
estimated  spectrum  at  these  frequencies  only. 

As  it  will  be  seen  in  the  sequel,  the  results  are 
satisfactory  even  in  the  presence  of  added  disturbances  in  the 
measurements.  The  recursive  use  of  blocks  of  data  has  the 
effect  of  smoothing  the  effect  of  disturbances  in  the 
estimates. 

D.   SIMULATION  RESULTS 

In  this  thesis,  we  try  inputs  with  different  frequencies 
to  test  the  identification  algorithm.  Also,  we  test  the 
behavior  when  the  input  is  the  combination  of  several 
different  frequencies  to  examine  the  effectiveness  of  the 
technique.  In  the  next  paragraph,  we  discuss  some  results 
obtained  by  running  the  program. 

1.   Original  System 

Figure  3 . 1  shows  the  magnitude  and  phase  of  frequency 
response  of  the  original  system.  It  shows  the  low  pass  nature 
of  the  system. 
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2.  Single  Frequency  Input  at  Arbitrary  Frequency  without 
Noise 

Figure  3 . 2  and  3 . 3  are  in  the  same  group  in  which  we 
tried  a  constant  input  with  magnitude  10.  As  predicted, 
figure  3 . 2  portrays  the  error  EPS  between  the  output  and 
predicted  output  which  decreases  and  approaches  zero. 
Furthermore,  in  figure  3.3,  the  frequency  response  of  the 
estimated  model  at  the  input  frequency  is  shown  to  be 
identical  to  the  one  of  the  original  system. 

In  a  second  set  of  experiments,  we  tried  several 
sinusoidal  inputs  with  different  frequencies.  First  of  all, 
we  used  an  input  signal  with  a  digital  frequency  tt/2  and 
obtain  the  results  shown  in  figures  3.4  and  3.5.  And  then,  we 
tried  another  input  with  frequency  tt/8  and  get  the  results 
shown  in  figures  3.6  and  3.7.  In  like  manner,  the  error 
approaches  to  zero  quickly  and  the  frequency  response  of  the 
estimated  model  at  the  input  frequency  is  the  same  as  the  one 
of  the  orignal  system. 

3.  Multiple  Frequency  Inputs  without  Noise 

In  a  third  set  of  experiments,  we  tried  the 
combination  of  four  frequency  inputs  at  arbitrary  frequencies 
to  obtain  the  results  shown  in  figures  3.8,  3.9,  and  3.10. 
The  output  in  the  time  domain  is  shown  in  figure  3.8.  In  this 
case,  the  error  between  true  and  predicted  outputs  approaches 


22 


zero  quickly  and  the   frequency  responses   at  the   input 
frequencies  are  still  identical. 

4.   Single  Frequency  Input  with  Noise 

In  this  section  we  show  results  of  similar  experiments 
with  random  measurement  noise.  The  noise  considered  is 
Gaussian  white  zero  mean,  i.e.,  independently  and  indentically 
distributed.  Its  standard  deviation  is  set  at  about  one  third 
of  the  output  magnitude. 

As  before  we  first  excite  the  system  with  a  constant 
input,  and  we  add  measurement  noise  to  the  output  of  the 
system.  The  results  are  shown  in  figures  3.11,  3.12,  and 
3.13.  From  figure  3.12,  even  though  we  find  that  the  error 
does  not  decay  to  zero  due  to  the  presence  of  noise,  the 
frequency  response  of  the  estimated  model  is  sufficiently 
close  to  the  one  of  the  original  system. 

In  a  second  set  of  experiments,  we  tried  sinusoidal 
inputs  again  with  noisy  measurements.  There  are  four  inputs 
with  different  digital  frequencies  (ll7r/16,7r/2 ,  7r/4,  and  tt/16) 
used  and  the  results  are  shown  in  fiqures  (3.14-3.22) 
respectively.  Then,  by  inspecting  the  figures,  the  error 
still  does  not  go  to  zero  due  to  the  presence  of  noise,  but 
the  estimated  frequency  response  is  sufficiently  close  to  the 
one  of  the  oriqinal  system. 
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5.   Multiple  Frequency  Inputs  with  Noise 

In  like  manner,  we  input  signals  with  multiple 
frequencies  to  test  the  alogrithm.  In  this  case,  we  used  four 
different  groups  of  inputs.  The  first  one  is  constant  and 
sinusoidal  (frequency  tt/2)  input.  The  2nd  one  is  a  constant 
and  sinusoidal  (frequency  n/4)  input.  The  3rd  and  4th  are 
sinusoidal  with  two  frequencies  (7r/2,7r/4)  and  (n/4  ,tt/8) 
respectively.  The  results  are  shown  in  figures  3.2  3  to  3.31. 
Then,  inspection  of  these  figures  reveals  that  the  error  in 
each  case  is  bounded  and  the  frequency  responses  of  the 
estimated  model  are  sufficiently  close  to  the  one  of  the 
original  system. 
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IV.   CONCLUSION 

In  this  thesis  we  introduced  an  approach  to  the 
identification  of  linear  models  based  on  the  frequency  domain 
formulation.  The  recursive  nature  and  the  use  of  FFT 
techniques  make  the  approach  attractive  for  on-line 
implementation . 

A  particular  aspect  we  addressed  is  the  identification  of 
the  frequency  response  of  the  system  by  approximating  the 
impulse  response  with  a  finite  sequence.  By  use  of  a 
simulated  example  we  have  explored  the  convergence  properties 
of  the  algorithm,  and  its  robustness  in  the  presence  of 
measurement  noise. 

Several  issues  remain  to  be  investigated,  mainly  the 
advantages  (if  any)  of  this  approach  in  comparison  with  more 
conventional  time  domain  estimation  techniques.  Although  this 
is  subject  of  future  research,  as  a  conjecture  we  can  say  that 
the  recursive  frequency  domain  approach  might  be  more  robust 
than  the  conventional  time  domain  in  the  cases  when  the  input 
signal  has  energy  concentrated  around  a  few  frequencies,  such 
as  the  case  of  periodic  excitation.  In  this  case  the 
processing  gain  of  the  DFT  should  give  a  better  performance  in 
the  presence  of  measurement  noise. 
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APPENDIX  A  :  MAINPROGRAM 


THESIS  MAIN  PROGRAM 


File  name 


THESIS4.M 


Advisor 
Student 
Date 


Roberto  Cristi 
Chao,  Chi-Shun 
Jan,  1,  1991 


%  Set  some  constant  number 


clg 

clear 

tfinal  =  80; 

dt  =  0.1; 

kmax  =  tfinal/dt; 

rand ( ' normal ' ) ; 

rand_mag  =  3  5; 

M  =  20; 

N  =  32; 

N2=  N/2; 


The  max.  time  to  run 

Sampling  time 

The  max.  number  of  signals 

Set  the  type  of  noise 

Set  the  gain  of  noise 

The  #  to  calculate  the  error 

The  #  of  signals  in  a  block 


%  Set  the 

yp(l)=0; 

yp(2)=0; 

y(l)=0; 

y(2)=0; 

u(l)=0; 

u ( 2 ) =0 ; 


initial  condition 


%  Produce  the  input  and  coresponding  output 
%  at  different  frequencies 
for  n=3:kmax 

%u(n)=  10; 

u(n)=  20*cos(pi*n*5/16) ; 
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%u(n)=  30*cos(pi*n/4) ; 
%u(n)=  40*cos(pi*n/8) ; 
%u(n)=  50*cos(pi*n/16) ; 
%u(n)=  60*cos(pi*n/32) ; 
%u(n)=  10+30*cos(pi*n/4) ; 
%u(n)=  40*cos(pi*n/8)+30*cos(pi*n/4) ; 

%u(n)=  10+2  0*cos(pi*n/2)+3  0*cos(pi*n/4)+4  0*cos(pi*n/8) 
%u(n)=  10+20*cos(pi*n/2)+30*cos(pi*n/8)  ... 
+40*cos (pi*n/8+50*cos (pi*n/16) ; 

yp(n)=0. 3*yp(n-l)+5. 0*u(n-l) ; 
y (n) =yp (n) +rand_mag*rand; 


end 


%  Produce  the  figures  of  input  and  output 

plot(u) ; 

title ( 'INPUT  U(n)  ' ) ; 

xlabel ( *  n ' )  ;   ylabel ( ' MAGNITUDE ' )  ; 

grid; 

pause; 

clg 

subplot(211) ,plot(yp) ; 

title ('YP(n)  ORIGINAL  OUTPUT  WITHOUT  NOISE  '); 

xlabel ( ■ n ' ) ;   ylabel ( ' MAGNITUDE ' ) ; 

grid; 

subplot (212) ,plot(y) ; 

title('Y(n)  OUTPUT  WITH  NOISE   u (n) =20cos (5*pi*n/16) 

y(n)=yp(n) ' ) ; 
xlabel ( ' n ' ) ;    ylabel ( ' MAGNITUDE ' ) ; 
grid; 

%meta  thesis44 
pause; 
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%  Produce  the  frequency  response  of  the  original  system 

%  by  using  the  "dbode"  function 

num=[5] ; 

den=[l,-0.3] ; 

w=linspace(0,pi, 17) ; 

[mag, phase] =dbode (num, den, w) ; 

clg 

subplot(211) ,plot(w,mag) ; 

title ( 'FREQUENCY  RESPONSE  OF  THE  ORIGINAL  SYSTEM  H(Z)'); 

xlabel( 'THETA' ) ; 

ylabel ( ' MAGNITUDE ' ) ; 

grid; 

subplot(212) , plot (w, phase)  ; 

title ( 'FREQUENCY  RESPONSE  OF  THE  ORIGINAL  SYSTEM  H(Z)1); 

xlabel( 'THETA' ) ; 

ylabel ( 'PHASE  (DEGREE)'); 

grid; 

pause; 


%   Create  the  matrix  H 
v= [ zeros ( 1 , N2 ) , ones ( 1 , N2 ) ] 
h=diag(v) ; 
F=fft(eye(N) ) ; 
H=F*h*inv(F) ; 


%  Set  the  zero  matrix 
THETAK  =  zeros(N,l); 
P0=10*eye(N) ; 
P=PO  ; 


59 


%  Use  the  FFT  and  RLS  (Recursive  Least  Square)  techniques 
%  to  estimate  frequency  response  of  the  model.  It  is  an 
%  interactive  and  iterative  loop 
for  k  =  2  :  M 

ykO  =  y(  (k-l)*N2+l  :  k*N2  )'; 

uk  =  u(  (k-2)*N2+l  :  k*N2  )'; 

yk= [ zeros (N2, 1) ;ykO] ; 

YK  =  fft(  yk  ) ; 

UK  =  fft(  uk  )  ; 

EPS (k) =0.0; 

for  i  =  1  :  N 

PHI=diag(H(i, : ) )*UK; 

[THETAK,P]=rls(THETAK,P,PHI,YK(i) ) ; 

EPS(k)=EPS(k)+(YK(i)-conj (PHI • ) *THETAK) ' *  ... 
(YK(i)-conj (PHI1 )*THETAK) ; 
end 

THETA(:,k)  =   THETAK  ; 

end 

EPS  =  sqrt(EPS)/N; 

%  Calculate  the  corresponding  horizontal  axis 
for  k=0:N/2 

wl(k+l)=2*pi*k/N; 
end 
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%  Produce  the  frequency  response  figures  of  estimated  model 

clg 

subplot (2 11) ,plot(wl,abs(THETA(l:17,M) ) ) ; 

title ( 'EST.  H(Z)  '  )  ; 

xlabel(»THETA     (  PI*K/(N/2)  )'); 

ylabel ( ' MAGNITUDE ' ) ; 

grid  ; 

subplot (2 12) , plot (wl, angle (THETA(1: 17, M) ) ) ; 

title ( 'EST.  H(Z)  ') ; 

Xlabel( 'THETA      (  PI*K/(N/2)  )'); 

ylabel ( 'PHASE  (RADIAN)'); 

grid; 

pause; 


%  Produce  the  error  figures 

clg 

displot(EPS) ; 

title ('EPS  (when  u(n) =20cos (5*pi*n/16) 

y (n)=yp(n)+35*rand   )•); 
xlabel( 'K') ; 
ylabel ( ' MAGNITUDE ' ) ; 
grid; 

Imeta  thesis44 
pause; 


%  Produce  the  frequency  response  of  the  estimated  model 

%  in  a  continuous  pattern. 

clg 

plot(w,mag,wl,abs(THETA(l:17,M) ) ) ; 

title('H(Z)  &  EST.  H(Z)  '); 

xlabel( 'THETA    (  PI*K/(N/2)  )'); 

ylabel ( ' MAGNITUDE '  )  ; 

grid; 

pause; 
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%  Store  the  value  of  the  response  of  the  estimated  model 
%  at  the  input  frequency. 

THETA_F  =  zeros(17,l); 


%THETA_F(1,1) 
%THETA_F(2,1) 
%THETA_F(3,1) 
%THETA_F(5,1) 
%THETA_F(9, 1) 
THETA_F (6,1) 


=  THETA(1,M) 
=  THETA(2,M) 
=  THETA(3,M) 
=  THETA(5,M) 
=  THETA(9,M) 
THETA(6,M) ; 


frequency  response  of  the  original  system  and 
model  at  the  input  frequency 


Produce  the 

the  response  of  estimated 

in  the  same  figure. 


clq 

plot (w,mag) ; 

title ('H(Z)  &  ESTIMATED  H(Z)  AT  THE  INPUT  FREQUENCY 

(WITH  NOISE) • ) ; 
xlabel( 'THETA    (  PI*K/(N/2)  )'); 
ylabel ( ' MAGNITUDE ' ) ; 
grid; 
hold  on 

displot (wl , abs (THETA_F) ) ; 
hold  off 
%meta  thesis44 
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APPENDIX  B  :  FUNCTION  TYPE  SUBROUTINE  RLS 


SUBROUTINE  TYPE  FUNCTION 


%  Advisor  :  Roberto  Cristi 


Student 
Date 


Chao,  Chi-Shun 
Jan,  1,  1991 


File  name 


RLS.M 


function  [x,p]=rls (x,p,phi,y) 


It  updates  the  estimate  x  and  its  covariance  p 
by  standard  recursive  least  squares. 


data: 


x  =  column  vector  (input,  output) 


p  =  square  matrix, 
y  =  scalar  (input) 
phi  =  column  vector 


positive  definite (input , output) 
(input) 


NOTE:   Due  to  numerical  approximations  it  might  be 

a  good  idea  to  check  the  matrix  p  for  positive 
def initeness. 


fact  =  1.0  +  phi'*p*phi; 

x  =  x  +  p*conj (phi) * (y-conj (phi' ) *x)/fact; 
p  =  p  -  p*conj (phi) *conj (phi' ) *p/ fact ; 
end  rls 
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